A quantum Monte Carlo study on the superconducting Kosterlitz-Thouless transition 
of the attractive Hubbard model on a triangular lattice 
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We study the superconducting Kosterlitz-Thouless transition of the attractive Hubbard model on 
a two-dimensional triangular lattice using auxiliary field quantum Monte Carlo method for system 
sizes up to 12 x 12 sites. Combining three methods to analyze the numerical data, we find, for the 
attractive interaction of U = — At, that the transition temperature stays almost constant within the 
band filling range of 1.0 < n < 1.4, while it is found to be much lower in the n < 1 region. 
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I. INTRODUCTION 

Discovery of superconductivity in layered materials 
or quasi-two-dimensional systems in the past several 
decades has brought up great interest in the physics 
of low dimensional superconductors. Among those are 
the high T c cupratesji a ruthenate Sr 2 RuC>4, 2 a cobal- 
tate Na^CoO^ • y H 2 Oj- MgB 2 ,A a heavy electron sys- 
tem CeCoIns^ and organic conductors such as (BEDT- 
TTF) 2 X (X is an anion) £ 

Theoretically, it is well known by Mermin- Wagner's 
theorem that no off-diagonal long range order takes place 
at finite temperature in purely two-dimensional(2D) 
systems^ In the case of pure 2D, superconducting tran- 
sition is expected to be of the Kosterlitz-Thouless (KT) 
type£ The superconducting KT transition has previ- 
ously been studied using finite temperature auxiliary field 
quantum Monte Carlo (AFQMC) technique for the at- 
tractive Hubbard model, that is, the Hubbard model with 
a negative on-site U, on a square lattice 9 - and on a tri- 
angular latticed Nowadays, a renewed interest for the 
KT transition on triangular lattices has arisen because 
unconventional superconductivity has been observed in 
materials having (anisotropic) triangular lattice structure 
such as Na x Co0 2 • yH 2 and (BEDT-TTF) 2 X. 

For the square lattice, it has been shown that the KT 
superconducting transition temperature T ~ O.lt, t is 
the hopping integral, for square lattice with n = 0.87 (n is 
the band filling) away from the half filling n = 1.0, where 
charge density wave (CDW) ordering takes place. 9 As for 
the triangular lattice, it has been concluded that Tq takes 
its maximum at around half-filling due to the absence 
of CDW ordering, while To is low near n — 1.4, where 
the density of states becomes large due to the van Hove 
singularity at n = 1.5. 

Since previous studies have been restricted to small 
system sizes, here we revisit the problem of the supercon- 
ducting KT transition of the attractive Hubbard model 
on a triangular lattice using AFQMC technique. Com- 
bining three methods of analyzing the numerical data, 
we find that the KT transition temperature stays almost 
constant within the band filling region of 1 < n < 1.4. Our 
results suggest that the estimation of the KT transition 
temperature using finite size scaling requires some cau- 
tion when the density of states near the Fermi level is 



small and the calculation is restricted to small system 
sizes. 



II. FORMULATION 
A. Model 

The Hamiltonian of the attractive Hubbard model is 
given in standard notation as 

i i,a 

where cJ CT (ci i(T ) is a fermion creation (annihilation) oper- 
ator at site i with spin a, and n,- )CT = c\ a Ci^ ai < i,j > de- 
notes a pair of nearest neighbors on square or triangular 
lattice having periodic boundary condition. The chemical 
potential fi controls the band filling n (average number of 
electrons per site). The hopping parameter t is taken as 
the unit of the energy and is set equal to unity through- 
out the study. As for the lattice structure, we mainly 
concentrate on the triangular lattice, but we also per- 
form calculation on the square lattice at the band filling 
of n = 0.87 in order to make comparison with previous 
studies and with the triangular lattice. The on-site at- 
traction U is fixed at U/t = —4.0 on both lattices, which 
also enables us to make comparison with the previous 
studies. 



B. Pairing Correlation Function 

Above the superconducting KT transition temperature 
Tq, the on-site (s-wave) pairing correlation function de- 
cays exponentially, 



< 4,T c U c Jd c id >~ ex P(- r /0- 



(2) 



By contrast, the pairing correlation function exhibits 
a power-law decay 



(3) 
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where rj = 0.25 for T -> T . In the actual AFQMC 
calculation, we calculate a summation of the pairing cor- 
relation function, P s , given as 

P s = (A^ + AA 1 ) (4) 



Due to this problem, here we identify To as the tem- 
perature at which P s for the largest two system sizes co- 
incide in the former method, while in the latter method, 
we search the values for To and A so that as many data 
points as possible fall on a single curve, although the 
data for small system sizes deviate from the curve at low 
temperatures. 



Due to the decaying behavior of the pairing correlation 
function mentioned above, P s increases as N is increased 
below To, while it remains constant above To. 

In the present study, finite temperature AFQMC is 
used to calculate P s m^i^ The number of Trotter de- 
composition slices L is chosen to satisfy the condition 
L > 12/3 (f3 = 1/T), so that At < 0.084/i is satisfied, 
where At = (3/L. 3000 - 30000 Monte Carlo sweeps, de- 
pending on the temperature and system size, have been 
taken to assure sufficiently small statistical errors. We 
have performed calculation on system sizes from 4 2 to 
12 2 sites. 

In order to obtain To from the AFQMC data of P s , 
we use three methods, two of which are based on finite 
size scaling, while the other one is a more straightforward 
method. 



C. Finite size scaling 

As shown in eq|3 the pairing correlation decays as 
r~ n below To for a large system size. Hence P s is pro- 
portional to Nx 2 ^ 71 with system size N = N 2 . In a finite 
size system, scaling variable N x /^ ,where £ is the corre- 
lation length, becomes more important. Taking this into 
account, the scaling hypothesis assumes the following be- 
havior for P s . 

P s (T,N x ) = N 2 x -if(N x /0 (6) 
^exp^AT-To) 1 / 2 ), (7) 

where A is a constant and f(x) is a certain scaling func- 
tion. Since £ — ► oo for T — * To, f(x) — /(0) regardless of 
the system size N at T = T . Thus, if we plot P s N- 2+r > 
as functions of T(or 0) for different system sizes, they 
should coincide regardless of the system size at T = T i& 

Another method to identify T is to plot P S N~ 2+V as a 
function of N x /£, where To and A are chosen so that all 
the data points fall on a single curve (namely P S N~ 2+11 = 
f(N x /£)) regardless of the system size£ 

However, as we shall see in the following, these scaling 
methods eqO turns out to suffer from finite size effects 
when the system size is small and/or the density of states 
at the Fermi level is small. In fact, in the former method 
mentioned above, P s for small system sizes do not coin- 
cide with those for large system sizes even at T = To. As 
for the latter method, P S N~ 2+V for small system sizes 
deviate from those for larger system sizes at low temper- 
atures. 



D. Extrapolation method 

Since the scaling method suffers from finite size ef- 
fects on some cases, we have adopted the third method, 
namely the extrapolation method. This method has 
been adopted for the study of the (real) superconducting 
transition temperature for the three dimensional attrac- 
tive Hubbard model, 13 Since P s should increase sharply 
at T = T) in the thermodynamic limit, it is reason- 
able to assume that To = liru/v— >oo , where Tj 
is the inflection point of P S (T) for the N site system. 

T (N) 

can be obtained by fitting P S (T) by an appro- 
priate function. In this study we use a gaussian form 
Am exp[— Bjy(T— Cn) 2 ] + Dn as a fitting function, where 
An-Dn are fitting parameters. Here, we find that 
scales as 1/iV" in the range of N considered in the present 

study. We may then obtain To by plotting as a func- 
tion of 1/N and linearly extrapolating it to the 1/N — > 
limit. 



III. RESULTS 
A. Finite size scaling 
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FIG. 1: P S N~ 1,75 plotted against (3 for the square lattice. 
Dashed lines are guides to the eyes. 

We first present finite size scaling analysis, where we 
identify the temperature at which TjiV" 1 - 75 coincides be- 
tween different system sizes. Before going into the re- 
sults for the triangular lattice, we show the result for the 
square lattice with n = 0.87, where a scaling analysis has 
been performed for system sizes up to 8 x 8 in ref|9|. In 
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a. n=0.8 



b. n=1.0 




n=1.2 



n=1.4 




FIG. 2: A plot similar to Fig^ but for triangu- 
lar lattice with several band filling. Band filling n is 
0.8(a), 1.0(b), 1.2(c), 1.4(d). 



FigE -RjiV" 175 is plotted as functions of (3 for system 
sizes from 4x4tol2xl2. We can see that the results 
for the largest two system sizes, 10 x 10 and 12 x 12 
merges at around (3 = 7. The results for smaller sys- 
tems do not cross or merge with each other within the 
present temperature range, but if we assume that the re- 
sults for small systems are strongly affected by finite size 
effects, we may adopt To/t ~ 1/7 ~ 0.13 for this band 
filling, which is roughly consistent with what has been 
concluded (T /t ~ 0.1) in ref® 

We now move on to the results for the triangular lat- 
tice. In Fig0 P s iV~ 1-7S is plotted as functions of [3. For 
n = 1.4, the results for the largest two systems, 10 x 10 
and 12 x 12 merges at around (3 = 6. The results for 
smaller systems do not cross or merge with each other 
within the present temperature range as in the case of 
the square lattice, but if we again assume that the re- 
sults for small systems are strongly affected by finite size 
effects, we may adopt T /t ~ 1/6 = 0.17 for this band 
filling. For n = 1.2, the results for the largest three sys- 
tem sizes cross at (3 ~ 5.3. Again assuming that results 
for smaller system sizes are affected by finite size effects, 
we may take Tb/t ~ 0.19 for this band filling. As for 
n = 1.0, although the results for 12 x 12 and 8x8 crosses 
at around (3 = 6, those of the largest two sizes do not 
cross with each other. Therefore, Tb/t cannot be esti- 
mated from these results at this band filling. The situa- 
tion is even worse in the case of n = 0.8, where P S N~ 1J5 
of the largest two system sizes do not even come close at 



low temperatures. Here again we cannot evaluate Tb/t 
for this band filling from these results. 



B. Extrapolation method 




FIG. 3: The AFQMC data of pairing correlation function P s 
for square lattice with n = 0.87. Dashed curves are the fitting 
results. 
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FIG. 4: plotted against 1/N for the square lattice with 

n = 0.87. The dashed line is a least squares fit. 



We now evaluate To using an alternative method, that 
is, the extrapolation method. Here again, we first show 
the result for the square lattice with n = 0.87. In Fig|3 
the raw P s data is plotted as functions of temperature 
for each system size. T 1 (Ar) extracted fr om these data are 
plotted against 1/N in Figg] The results for N = 4 2 
turn out to be too much affected by finite size effects, so 
we have omitted these data in the extrapolating process. 
T obtained by extrapolating the results to 1/N — > is 
To/t ~ 0.15, which is a little bit higher than, but in fair 
agreement with the value obtained by the scaling method. 

We move on to the triangular lattice. In Fig0 the 
raw P s data as well as the fitting curves are plotted as 
functions of T for each system size and band filling. (We 

do not show the fitting curves for n = 0.8 since 

turns out to be negative.) obtained fromthese data 

are plotted against 1/N for each band filling in FigO To 
obtained by extrapolation is 0.17, 0.18, 0.17 for n = 1.4, 
n = 1.2, n = 1.0, respectively. For n = 1.4 and n = 1.2, 
Tq obtained using the scaling method coincides fairly well 
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c. n=1.2 



d. n=1.4 




FIG. 5: A plot similar to Fig|^]but for the triangular lattice 
with n = 0.8(a), 1.0(b), 1.2(c), 1.4(d). Dashed lines are fitting 
results. 
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FIG. 6: A plot similar to Fig[I]but for the triangular lattice 
with n = 1.0, 1.2, 1.4. 
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FIG. 7: The band filling dependence of the KT transition 
temperature To. Solid and dashed lines are guides to the 
eyes. 
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FIG. 8: PsN' 1 - 75 plotted against N x /£ for the triangular 
lattice. 



with the values obtained here. For n — 0.8, P s turns out 
to decrease upon increasing the system size from 8 x 8 to 
12 x 12, which seems to indicate that there is no symptom 
of KT transition within the temperature range (> O.li) 
investigated in the present study, so To, if any, should 
be lower than O.li. In Fig0 T obtained by the present 
method is plotted against the band filling. 



C. Justification by scaling 

The KT transition temperature determined by the 
above method can be further justified by checking 
whether P s actually scales as P s = N^ r5 f(N x /£) with 



£ = exp(A/(T - To) 1 / 2 )). In FigEl we plot P.JV" 1 - 78 as 
functions of N x /£ for all the system sizes. We find for 
n = 1.0, 1.2, and 1.4 that the numerical results for 8x8, 
10 x 10 and 12 x 12 fall on a single curve by adopting 
the To obtained above and choosing appropriate values of 
A. The results for smaller systems also fall on the same 
curves at high temperatures (namely for large N x /£), but 
at low temperatures they start to deviate due to finite 
size effects. As for n = 0.8, if A = 0.3 and T = 0.05 
are chosen, the results for 12 x 12 and 10 x 10 fall on a 
similar curve. These results further justify the values of 
To obtained in the preceding sections. Similar results are 
also obtained for the square lattice as seen in Fig|5] 
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FIG. 9: A plot similar to Fig|H|for the square lattice with 
n = 0.87. 



IV. DISCUSSION 

A. Comparison of the methods for determining To 

Our results show that the scaling method using the 
crossing point of PgN' 1 - 75 is strongly affected by finite 
size effects, and that it is difficult to know from the begin- 
ning the system size necessary to obtain an accurate To. 
The present analysis suggests that this method seems to 
work well at band fillings where the density of states at 
Ep is relatively large, namely at n = 1.2 and n = 1.4 in 
the present case. There, the agreement with the results 
of the extrapolation method is also good. This may be 
because the discreteness of the energy levels due to the 
finite system size is small when the density of states is 
large. 



B. Correlation between To and the density of states 

In order to look into a possible correlation between To 
and the density of states at the Fermi level, we compare 
the present results of To with the transition tempera- 
ture obtained by mean field approximation (Figflip. Al- 
though the values of T^ F themselves are much larger 
than T), wc find a similar band filling dependence, 
namely, T™ F is almost constant for 1 < n < 1.5, and 
smaller for n < 1. The origin of this similarity is not 
clear, but this does seem to suggest that T is roughly 
correlated with the density of states near the Fermi level 
at least for the value of U (= —At) adopted in the present 
study. It would be an interesting future study to ana- 
lyze this point from the viewpoint of the study by Timm 
et al, where the relation between the mean field T c and 
the KT superconducting transition temperature in the 
repulsive Hubbard model has been investigated by cal- 
culating the superfluid density as a function of temper- 
ature and combining it with the Berezinskii-Kosterlitz- 
Thouless theory>ii 

T) for 1 < n < 1.4 is still larger than Tq for the square 
lattice at n = 0.87 despite the fact that the density of 




band filling n 

FIG. 10: The density of states of the square and the trian- 
gular lattices on the tight binding model. 
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FIG. 11: The band filling dependence of the T C MF for the 
triangular lattice obtained by applying BCS mean field ap- 
proximation. 



states near the Fermi level for the square lattice with 
n = 0.87 is larger than for the triangular lattice with 
n = 1.0 (see FigUUJl. 

If T) is indeed positively correlated with the density 
of states near the Fermi level, this "inversion" of To be- 
tween the square and the triangular lattice may be due to 
the absence of charge density wave ordering in the latter 
lattice due to frustration, as discussed previously^ 



V. CONCLUSIONS 

In the present study, we have investigated the su- 
perconducting KT transition of the attractive Hubbard 
model on a two-dimensional triangular lattice using aux- 
iliary field quantum Monte Carlo method for system sizes 
up to 12 x 12 sites. Combining three methods to analyze 
the numerical data for the pairing correlation function, 
we find that the transition temperature stays almost con- 
stant within the band filling range of 1.0 < n < 1.4, while 
it is found to be much lower in the n < 1 region. Among 
the three methods, the extrapolation method is found to 
work well regardless of the band filling, while the methods 
relying on finite size scaling require some caution when 
the density of states near the Fermi level is small and the 
calculation is restricted to small system sizes. 
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